Dynamical compressibility of dense granular shear flows 
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It has been conjectured by Bagnold [lj that an assembly of hard non-deformable spheres could 
behave as a compressible medium when slowly sheared, as the average density of such a system 
effectively depends on the confining pressure. Here we use discrete element simulations to show 
the existence of transverse and sagittal waves associated to this dynamical compressibility. For this 
purpose, we study the resonance of these waves in a linear Couette cell and compare the results 
with those predicted from a continuum local constitutive relation. 



Acoustics in granular media is a rapidly developing 
field that is attractive both from the fundamental point 
of view and for its applications in material science, soil 
mechanics and geophysics It provides a unique tool 

to probe the mechanical response of such materials, to 
detect heterogeneities, aging or avalanche precursors. In 
the static case, grain elasticity is the only restoring force; 
elastic waves present most of the known behaviors of me- 
chanical waves: velocity dispersion, non-linearity associ- 
ated to contact geometry, scattering associated to disor- 
der, anisotropy, aging, mode coupling, wave-guiding in 
thepresence of macroscopic heterogeneities, fragility, etc 
Recent studies 0-11 1 have focused on the failure of 
the mean field description of elasticity when approaching 
the jamming transition separating the solid- like from the 
liquid-like regime. The importance of non-affinc displace- 
ments in this limit was revealed by the structure of the 
vibration spectrum, which exhibits an anomalous excess 
of low frequency 'soft' modes fl2l - [l4l ]. In the following 
we wish instead to investigate wave propagation on the 
liquid side of the jamming transition. 

Wave propagation in granular shear flows have never 
been studied for itself. It has been revealed indirectly by 
instabilities, in situations where acoustic waves are spon- 
taneously emitted during a dense granular flow (see Q 
for a review): when sufficiently dry, avalanches flowing 
at the surface of a du ne p roduce a loud sound known 
as the 'song of dunes' [15l4l8l | : vibrations are also pro- 
duced during the discharge of a granular flow from a silo 
or any elongated tube [19(; soft mono-disperse particles 
flowing on an inclined plane near the angle of repose ex- 
hibit spontaneous oscillations at the resonant frequency 
of elastic waves 2fj, 21 1 • These phenomena have mostly 
been associated to the elastic deformations of the grains. 
However, Bagnold [l[ and followers M HI have hy- 
pothesized the existence of vibrations driven by the cou- 
pling between normal stress and shear rate. 

In this Letter we theoretically demonstrate the exis- 
tence of such waves in dense slow shear flows of hard 
non-deformable grains. We show that these waves are 
due to a compressibility of dynamical origin, whose scal- 
ing with the volume fraction <j> is related to the non-affinc 
cooperative motion of the grains. The dispersion relation 
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FIG. 1: (a) Schematic of the numerical set-up. (b) Typical 
velocity signal as a function of time t, in a simple shear flow 
at (f> ~ 0.78, H/d ~ 40, and yd ~ 0.1-y/V/p. To improve the 
signal-to-noise ratio, we have computed the projection u z of 
the vertical velocity field on the mode sm(nz/ H). (c) Auto- 
correlation functions C(jr) of the signal u z (t), for different 
values of the parameters K n , n t and e. The first maximum 
gives the spontaneous frequency uj* /2n. 

presents three branches, two of which becoming non- 
dissipative when approaching the jamming point. Our 
analytical predictions are confirmed by numerical sim- 
ulations performed with quasi-rigid grains, which show 
no dependence on their elasticity nor on their restitution 
coefficient. 

Numerical model - We consider a two-dimensional sys- 
tem constituted of ~ 2 • 10 3 circular particles of mass rrii 
and diameter di, with a ±50 % polydispersity (Fig. QJ,). 
The shear cell is composed of two rough walls distant 
by H and moving along the x-dircction at a constant 
velocity. Simulations are therefore performed at a con- 
stant, imposed volume fraction 0. Periodic boundary 
conditions are applied along x. The particle dynamics is 
integrated using the Verlet algorithm. The particles are 
submitted to contact forces modeled as viscoelastic forces 
along the normal direction an d as a Coulomb friction 
along the tangential direction 



22j-[24|. Quantities used 



in the model are expressed in terms of the grain density 
p, of the average pressure V , and of the mean grain diam- 
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FIG. 2: (a) Spontaneous angular frequency uu* as a function 
of (4> c — (f>). Solid line: prediction for the resonant frequency 
2nc/H deduced from Eq. [7] Insets: Black solid line: power 
spectral density D(lj) of the projected vertical velocity u z . 
Gray dotted line: best fit by a gaussian. (b) Scaling of the 
pressure with the volume fraction. Solid line: best fit by the 
form f(4>) given by Eq. ||2j. 

eter d. In this system of units, the normal spring constant 
K n is chosen sufficiently large (between 10 3 and 10 4 ) to 
reach the rigid asymptotic regime. The Coulomb friction 
coefficient is chosen equal to fi p = 0.4 and the tangential 
spring constant Kt = 0.5k„. The influence of viscoelas- 
tic parameters is encoded in the restitution coefficient e, 
which is varied between 0.1 and 0.9. The average shear 
rate 7 is kept constant during all simulations. 

Spontaneous and forced oscillations - Looking at the 
time fluctuations of the velocity field (Fig. QTj), sim- 
ple shear flows are found to present spontaneous vibra- 
tions. As seen from the velocity autocorrelation function 
(Fig. [TJ:) and the power spectral density (Fig. [2|) of the 
velocity signal, one observes that this vibration is not a 
broadband noise but oscillations at a well defined angular 
frequency u>*. Keeping the shear rate 7 constant and ap- 
proaching the jamming point <p c , tu* is found to diverge as 
(4> c — 0)~ 3 / 2 (Fig. [2^). The coherence time can be deter- 
mined from the spectrum, which presents a narrow peak 
whose width at half maximum is around w*/4. (Fig. [2]). 
Furthermore, as (f> — > <p c , it develops a tail towards small 
frequencies, which is reminiscent of the spectral density 
of elastic modes on the solid side of jamming, for cf> > cf> c . 

To analyze the vibration mode by mode, we have 
performed simulations in which each grain is submitted 
to a force along z of the form f z sin(wt) sin(fcz), with 
k = tt/H. We analyze the response of the system to this 
forcing with a simple 'lock-in amplifier': we compute the 
Fourier component of the velocity field at the angular 
frequency u) and the wavenumber k. The amplitude of 
the forcing is chosen small enough to remain in the linear 
regime. The velocity response is then proportional to the 



forcing f z and presents a resonance, characterized by a 
maximum of amplitude and by a vanishing phase (Fig. [3]). 
The resonant angular frequency coincides with the angu- 
lar frequency uj* of spontaneous oscillations (Fig. [TJ). 

The resonant curves (Fig. [3J and the autocorrelation 
function of spontaneous vibrations (Fig. [T};) turn out to 
be insensitive to the values of the restitution coefficient e 
and of the spring constants. The lack of dependence on 
e shows that the dense flow docs not behave like a gran- 



ular gas [251 - 4281 ] . whose compressibility is controlled by 
binary collisions. The nature of these oscillations is also 
fundamentally different from those exhibited in [20|, [2l| 
which are specifically observed for mono-disperse grains 
and whose frequency scales The vibration modes 

of a system of hard sheared particles are therefore differ- 
ent from both elastic modes in solids and sound waves in 
granular gases. 

Continuum model - To shed light on the origin of 
these spontaneous oscillations, we consider a descrip- 
tion of granular material as a continuum compressible 
medium characterized by the volume fraction <p and the 
velocity field u. Under the assumption that the rheology 

H 



is local 



301 , dimensional analysis constrains the form 
r e relation betweer 
strain rate 7^ = (diiij + djUi) to: 



of the constitutive relation between the stress Oij and the 
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where fj,(<f>) is the effective friction coefficient, /(</>) the 
rescaled dilatancy, and g(<p) is the rescaled bulk viscosity. 
As shown in Fig. [2Jd, f(<f>) diverges as 

m = Ti (2) 



at the critical volume fraction <p c 31 - 34| ■ The friction is 
regular in the limit <p — > (f> c and is well approximated by 



= Mc 



Sfi 



l + ab/{(/) c - 4>)' 



(3) 



We hypothesize that the rescaled bulk viscosity is reg ular 
as well: g — > g c as <p (f) c - In a previous study [34|, we 
have obtained <j> c ~ 0.817, a — 0.36, b = 0.33 and S/j, = 
0.57 in the 2D case, using homogeneous shear flows, g has 
never been measured so far. Importantly, the equations 
are the same in both 2D and 3D. 

The compressible Navicr-Stokcs equations are given by 



Du _ 

) = Ver + F and 

Dt 



^ = -^(V-u), (4) 



where F is a bulk force. The base state is a homogeneous 
steady-state shear flow, and from now on we denote its 
volume fraction by (f>, its shear rate by 7 and its pressure 
by V = /(0)p<i 2 |7| 2 . Linear perturbations in volume 
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FIG. 3: (Color online) Amplitude and phase of the velocity response along the z axis (a) and the x axis (b). Solutions of 
Eqs. ((9J, with forcing, are displayed in solid (with a bulk viscosity g = 2) and dotted lines (for g = 0). The velocity field is 
projected onto the form Uj sin(fcz) sin(cj£ + aj), where «j is the amplitude and Oi the phase shift with respect to the forcing. 
The various symbols correspond to different restitution coefficients, different spring constants and different amplitudes of the 
forcing. The simulation conditions are the same as in Fig. [TJ with FJ ef ~ [j?) fz Ct — 0.025 V/d. 



fraction and velocities around the base state will be de- 
noted by 1 and u\. Plane waves can only propagate 
along the y and z axis as a plane wave emitted along 
x would rotate due to convection by the base velocity 
field u x = jz. This effect can be alternatively seen as a 
refraction by an index gradient. 

Let us consider first the case of waves propagating 
along the transverse axis (in 3D), of the form exp(iujt + 
iky). In order to induce compression, the wave must have 
a velocity component Uy along y, which corresponds to a 
strain perturbation 7^ = 2ikUy. The equation of motion 
projected along y and the continuity equation then read 



The imaginary part of the dispersion relation gives the 
quality factor Q which, at the jamming point scales as: 



Q = - 
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(8) 



This term diverges as {4> c — 0) _1 , which means that the 
propagation of this transverse mode becomes asymptot- 
ically non-dissipativc. 

We now consider modes propagating along z, of the 
form cxp(iujt + ikz). The perturbation of the strain 
rate components are jl z = iku\ and ^\ z = 2iku\. The 
Navier-Stokes equations then read 



-cj)ikv} y , 
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In this equation and in the followings, the functions /, 
/', g and [i are evaluated at the volume fraction <j) of the 
base state. One finally obtains the dispersion relation 
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When the imaginary term can be neglected in front of the 
real one, the propagation is found non-dispersive, with a 
speed 



This system can be easily solved to derive the disper- 
sion relation, which presents two branches. Approaching 
the jamming transition, as (f> — > <p c , the first of these 
two modes propagates at a speed equal to c (Eq. [7|); it 
presents an elliptical polarization described by 
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At a fixed shear rate, the speed of sound is therefore 
found to diverge at the jamming point as (<p c — 0)" 3 / 2 . 



u z U 

and a quality factor diverging at the jamming point as 
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The second mode is critically damped and presents a dif- 
fusive dispersion relation (k 2 oc ui) close to jamming: 
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It becomes linearly polarized along a; as 
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As it is highly dissipative, this shear mode cannot be 
observed. 

Fig. [2ji shows that the spontaneous angular frequency 
coincides, without any adjustable parameter, to the res- 
onant mode 2irc/H deduced from Eq. ([7|). Adding a si- 
nusoidal bulk force on the right hand side of Eqs. (J9j) 
and applying a no-slip boundary condition at the walls, 
one obtains the predictions for the resonance curves, with 
a single adjustable parameter: the second viscosity. As 
shown in Fig. [31 the resonant frequency, the phase curve 
and the tails of the amplitude response are largely in- 
sensitive to the value of g. The fit of the amplitude 
around the resonance for different volume fractions gives 
g(<j>) ~ g c ~ 2. 

Concluding remarks - Both transverse and sagittal 
modes are directly controlled by the dynamical compress- 
ibility, which results from the non-affine particle motion 
during shearing. Close to the jamming point, the sta- 
tistical properties of grain trajectories are remarkably 
insensitive to the nature of the dissipative mechanisms 
31, [32! and to the dynamical regime considered (iner- 
tial or overdamped) 34|. As steric effects dominate, 
cooperative non-affine motions mostly depend on (f>. To 
move by an average distance d along x, a grain actually 
makes an erratic motion whose average length diverges 



as 



d(cf> c — <j>) [3l|,l35|. By definition of the shear rate, 
such a displacement occurs at a rate ~ 7. The dynamical 
pressure, controlled by inertia, scales as the squared ve- 
locity fluctuations and therefore as V ~ pj 2 d 2 ((j) c — (f))~ 2 
(see Fig. [2Jd). Let us emphasize that this scaling law is 
different from that obtained in the unsheared case, ther- 
malized j3(| or not [37j |. for which V ~ (0 C — t/))" 1 and 
therefore c ~ (<f) c — (f))^ 1 ■ 

The propagative modes exhibited in this Letter are 
therefore specific of the sheared dense liquid granular 
regime: they neither depend on the restitution coef- 
ficient nor on the grain elasticity and the scaling law 
obeyed by the speed of sound c is a direct consequence 
of cooperative effects close to jamming. As these modes 
become asymptotically non-dissipative at the jamming 
point, they should be observed experimentally and could 
be used to probe the bulk of granular flows. 
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